home *** CD-ROM | disk | FTP | other *** search
/ Workbench Add-On / Workbench Add-On - Volume 1.iso / BBS-Archive / Dev / ace23.lha / PRGS.lha / Astronomy / Jovian.b < prev    next >
Text File  |  1994-10-02  |  12KB  |  541 lines
  1. { Displays the positions of the Galilean
  2.   satellites relative to Jupiter for a 
  3.   given date and period of time.
  4.  
  5.   The view is as it would appear through
  6.   an inverting telescope in the Southern
  7.   Hemisphere.
  8.  
  9.   All dates and times must be entered in 
  10.   Universal (Greenwich Mean) Time.
  11.  
  12.   Method taken from Jean Meeus' 
  13.   "Astronomical Algorithms", ch 43.
  14.  
  15.   Author: David J Benn
  16.     Date: 11th,12th,13th,17th,26th July 1993,
  17.       18th December 1993,
  18.       7th March 1994,
  19.       14th August 1994 }
  20.  
  21. CONST radconv=57.295779
  22. CONST xorigin=320!, yorigin=116!
  23. CONST radius=5
  24. CONST true=-1&, false=0&
  25. CONST JovianColor=2
  26. CONST black=0
  27. CONST xscale=10, yscale=7.5
  28.  
  29. SINGLE   JDE
  30. SINGLE   d,V,M,N,J,A,B,K,R,rr,delta,psi,DE
  31. SINGLE   first_u1,first_u2,first_u3,first_u4
  32. SINGLE   u1,u2,u3,u4
  33. SINGLE   G,H
  34. SINGLE   r1,r2,r3,r4
  35. LONGINT  first
  36. SHORTINT moon
  37.  
  38. DIM x(4),y(4),lastx(4)
  39.  
  40. SUB SINGLE decimal_hours(hrs,mins,secs)
  41. '..return decimal hours
  42.  
  43.   decimal_hours = hrs + mins/60 + secs/3600
  44. END SUB
  45.  
  46. SUB SINGLE JulianDay(dy,mn,yr) 
  47. SINGLE m1,y1,a,b,c,d,dj
  48. '..This routine calculates the number of days elapsed 
  49. '..since the epoch 1900 January 0.5 (ie: 1200 GMT, 31st Dec 1899). 
  50.  
  51.  if yr = 0 then 
  52.    dj=-1  '..error 
  53.  else 
  54.    m1=mn : y1=yr : b=0 
  55.  
  56.   if y1 < 1 then ++y1 
  57.   if mn < 3 then m1=mn+12 : --y1 
  58.  
  59.   if y1 > 1582 or mn > 10 or dy >= 15 then  
  60.     a=int(y1/100) : b=2-a+int(a/4) 
  61.     c=int(365.25*y1)-694025 
  62.     if y1 < 0 then c=fix((365.25*y1)-0.75)-694025 
  63.     d=int(30.6001*(m1+1))
  64.     dj=b+c+d+dy-0.5 
  65.   else 
  66.     if (y1<1582 or (y1=1582 and mn<10) or (y1=1582 and mn=10 and dy<5)) then
  67.       c=int(365.25*y1)-694025 
  68.       if y1 < 0 then c=fix((365.25*y1)-0.75)-694025 
  69.       d=int(30.6001*(m1+1)); dj=b+c+d+dy-0.5 
  70.     else 
  71.       dj=-1  '..error
  72.     end if
  73.   end if 
  74.  end if
  75.  
  76.  JulianDay = dj  '..Return Julian Date (error = -1)
  77. END SUB
  78.  
  79. SUB STRING time_from_day_fraction(SINGLE fd)
  80. '..return time string from day fraction.
  81. SINGLE hrs,mins
  82.   
  83.   hrs = 24*fd
  84.   mins = 60*(hrs-fix(hrs))
  85.   hr$ = str$(fix(hrs))
  86.   min$ = str$(fix(mins))
  87.   hr$ = right$(hr$,len(hr$)-1)
  88.   min$ = right$(min$,len(min$)-1)
  89.   if len(hr$)=1 then hr$="0"+hr$
  90.   if len(min$)=1 then min$="0"+min$
  91.  
  92.   time_from_day_fraction = hr$+":"+min$
  93. END SUB
  94.  
  95. SUB STRING date_and_time(dj!)
  96. '..This routine converts the number of (Julian) days since 
  97. '..1900 January 0.5 into the calendar date and time.
  98. SINGLE a,b,c,d,g,i,fd
  99.  
  100.  d=dj!+0.5 : i=int(d) : fd=d-i 
  101.  
  102.  if fd = 1 then fd=0 : ++i 
  103.  
  104.  if i > -115860 then 
  105.    a=int((i/36524.25)+9.9835726e-1)+14 
  106.    i=i+1+a-int(a/4) 
  107.  end if 
  108.  
  109.  b=int((i/365.25)+8.02601e-1) 
  110.  c=i-int((365.25*b)+7.50001e-1)+416 
  111.  g=int(c/30.6001) : mn=g-1 
  112.  dy=c-int(30.6001*g)+fd : yr=b+1899 
  113.  if g > 13.5 then mn=g-13 
  114.  if mn < 2.5 then yr=b+1900 
  115.  if yr < 1 then --yr 
  116.  
  117.  '..return a date string (whole days only) 
  118.  dy$=str$(int(dy)) : if dy < 10 then dy$="0"+right$(dy$,1)
  119.  dy$=right$(dy$,2)
  120.  mn$=str$(int(mn)) : if mn < 10 then mn$="0"+right$(mn$,1)
  121.  mn$=right$(mn$,2)
  122.  yr$=str$(int(yr))
  123.  yr$=right$(yr$,4)
  124.  
  125.  date_and_time = dy$+"-"+mn$+"-"+yr$+"   "+time_from_day_fraction(fd)
  126. END SUB 
  127.  
  128. SUB SINGLE in_range(n)
  129. '..ensure n is in the range 0..360 degrees
  130.   if n<0! then 
  131.     in_range = 360! + (n mod 360!)
  132.   else
  133.     if n>360! then in_range = n mod 360!
  134.   end if
  135. END SUB
  136.  
  137. SUB SINGLE sinh(x)
  138. '..return hyperbolic sine of x
  139.   sinh = (exp(x)-exp(-x))/2!
  140. END SUB
  141.  
  142. SUB JovianEphemeris(SINGLE JDE)
  143. SHARED d,V,M,N,J,A,B,K,R,rr,delta,psi,DE
  144. '..calculate circumstances of Jupiter at JDE
  145.  
  146.   '..days since 2000 January 1, 12h
  147.   d = JDE - 36525.0
  148.  
  149.   '..argument for long-period term in motion of Jupiter    
  150.   V = in_range(172.74 + 0.00111588*d)    
  151.  
  152.   '..mean anomalies of Earth and Jupiter
  153.   M = in_range(357.529 + 0.9856003*d)
  154.   N = 20.02 + 0.0830853*d + 0.329*sin(V/radconv)
  155.  
  156.   '..difference between mean heliocentric 
  157.   '..longitudes of Earth and Jupiter
  158.   J = in_range(66.115 + 0.9025179*d - 0.329*sin(V/radconv))
  159.  
  160.   '..equations of the center of Earth and Jupiter
  161.   A = 1.915*sin(M/radconv) + 0.02*sin((2*M)/radconv)
  162.   B = 5.555*sin(N/radconv) + 0.168*sin((2*N)/radconv)
  163.   K=J+A-B
  164.  
  165.   '..radius vector of Earth
  166.   R = 1.00014 - 0.01671*cos(M/radconv) - 0.00014*cos((2*M)/radconv)
  167.  
  168.   '..radius vector of Jupiter
  169.   rr = 5.20872 - 0.25208*cos(N/radconv) - 0.00611*cos((2*N)/radconv)  
  170.  
  171.   '..Earth-Jupiter distance
  172.   delta = ABS(SQR(rr*rr + R*R - 2*rr*R*cos(K/radconv)))
  173.  
  174.   '..phase angle (Earth-Jupiter-Sun)
  175.   sin_of_psi = (R/delta)*sin(K/radconv)
  176.   psi = sinh(sin_of_psi)*radconv
  177.  
  178.   '..longitudes of central meridian in systems I and II
  179.   w1 = in_range(210.98 + 877.8169088*(d-(delta/173!)) + psi - B)
  180.   w2 = in_range(187.23 + 870.1869088*(d-(delta/173!)) + psi - B)
  181.  
  182.   '..heliocentric longitude
  183.   lambda = 34.35 + 0.083091*d + 0.329*sin((V+B)/radconv)
  184.  
  185.   '..planetocentric declination
  186.   DS = 3.12*sin((lambda+42.8)/radconv)
  187.   DE = DS - 2.22*sin(psi/radconv)*cos((lambda+22!)/radconv)
  188.   DE = DE - 1.3*((rr-delta)/delta)*sin((lambda-100.5)/radconv)  
  189. END SUB
  190.  
  191. SUB AngleFromInfConj
  192. '..calculate angle from inferior conjunction
  193. SHARED d,delta,psi,B
  194. SHARED first_u1,first_u2,first_u3,first_u4
  195. SHARED u1,u2,u3,u4
  196. SHARED G,H
  197.  
  198.   deltaterm = d - (delta/173)
  199.  
  200.   first_u1 = in_range(163.8067 + 203.4058643*deltaterm + psi - B)
  201.   first_u2 = in_range(358.4108 + 101.2916334*deltaterm + psi - B)
  202.   first_u3 = in_range(5.7129 + 50.2345179*deltaterm + psi - B)
  203.   first_u4 = in_range(224.8151 + 21.4879801*deltaterm + psi - B)
  204.  
  205.   '..correct for mutual perturbations
  206.   G = 331.18 + 50.310482*deltaterm
  207.   H = 87.4 + 21.569231*deltaterm
  208.  
  209.   u1 = first_u1 + 0.473*sin((2*(first_u1-first_u2))/radconv)
  210.   u2 = first_u2 + 1.065*sin((2*(first_u2-first_u3))/radconv)
  211.   u3 = first_u3 + 0.165*sin(G/radconv)
  212.   u4 = first_u4 + 0.841*sin(H/radconv)
  213. END SUB
  214.  
  215. SUB DistFromCenter
  216. '..calculate distance of each satellite from
  217. '..center of Jupiter
  218. SHARED first_u1,first_u2,first_u3,first_u4
  219. SHARED r1,r2,r3,r4
  220. SHARED G,H
  221.   
  222.   r1 = 5.9073 - 0.0244*cos((2*(first_u1-first_u2))/radconv)
  223.   r2 = 9.3991 - 0.0882*cos((2*(first_u2-first_u3))/radconv) 
  224.   r3 = 14.9924 - 0.0216*cos(G/radconv)
  225.   r4 = 26.3699 - 0.1935*cos(H/radconv)
  226. END SUB
  227.  
  228. SUB calc_x_y(SHORTINT n)
  229. '..calculate rectangular coordinates 
  230. '..of four satellites
  231. SHARED x,y
  232. SHARED r1,r2,r3,r4
  233. SHARED u1,u2,u3,u4
  234. SHARED DE
  235.  
  236.   case
  237.     n=1 : r=r1:u=u1
  238.     n=2 : r=r2:u=u2
  239.     n=3 : r=r3:u=u3
  240.     n=4 : r=r4:u=u4
  241.   end case 
  242.  
  243.   x(n) = r*sin(u/radconv)
  244.   y(n) = -r*cos(u/radconv)*sin(DE/radconv)   
  245. END SUB
  246.  
  247. SUB LONGINT moving_east(SHORTINT moon)
  248. '..return true if moon is moving east.
  249. SHARED lastx,x  
  250.  
  251.   if lastx(moon) > x(moon) then 
  252.     moving_east = true
  253.   else
  254.     moving_east = false
  255.   end if
  256. END SUB
  257.  
  258. SUB LONGINT in_disk_region(SHORTINT x,SHORTINT y)
  259. '..return true or false according to whether
  260. '..a satellite is in the region of the disk
  261. '..of Jupiter.
  262.  
  263.   if point(x-1,y)=JovianColor and point(x+1,y)=JovianColor then
  264.     in_disk_region = true
  265.   else
  266.     in_disk_region = false
  267.   end if    
  268. END SUB
  269.  
  270. SUB LONGINT behind_disk(SHORTINT xcoord,SHORTINT ycoord,SHORTINT moon)
  271. '..return true or false according to whether
  272. '..a satellite is behind the disk of Jupiter.
  273.  
  274.   if in_disk_region(xcoord,ycoord) and moving_east(moon) then
  275.     behind_disk = true
  276.   else
  277.     behind_disk = false
  278.   end if   
  279. END SUB
  280.  
  281. SUB RemoveMoons
  282. SHARED x,y
  283. SHORTINT moon,xx,yy,colr
  284.  
  285.   '..clear moons
  286.   for moon=1 to 4
  287.     xx = xorigin+x(moon)*xscale
  288.     yy = yorigin-y(moon)*yscale
  289.     
  290.     if in_disk_region(xx,yy) then
  291.       '..in transit across disk or behind it
  292.       colr = JovianColor
  293.     else
  294.       colr = black
  295.     end if
  296.     pset (xx,yy),colr
  297.   next
  298. END SUB
  299.  
  300. SUB ShowJovianSpace
  301. '..display Jupiter and the Galilean satellites
  302. SHARED x,y
  303. SHORTINT xx,yy,moon
  304.  
  305.   '..plot moons
  306.   for moon=1 to 4
  307.     xx = xorigin+x(moon)*xscale
  308.     yy = yorigin-y(moon)*yscale
  309.     if not behind_disk(xx,yy,moon) then pset (xx,yy),moon
  310.   next
  311.  
  312.   '..draw Jupiter
  313.   circle (xorigin,yorigin),radius,JovianColor
  314.   paint (xorigin,yorigin),JovianColor
  315. END SUB
  316.  
  317. SUB DisplayData
  318. '..display Jupiter/satellite data 
  319. CONST startcol=15
  320. SHARED JDE,x,y,rr,delta
  321. SHORTINT moon  
  322.  
  323.   '..Date & Time
  324.   locate 3,10
  325.   color 1,0
  326.   print date_and_time(JDE);" UT"
  327.   
  328.   '..Earth-Jupiter distance  
  329.   locate 4,10
  330.   color 4,0
  331.   print "Earth-Jupiter Distance (AU): ";
  332.   color 5,0
  333.   print delta;"    " 
  334.  
  335.   '..Jupiter's Radius Vector
  336.   locate 5,10
  337.   color 6,0
  338.   print "  Sun-Jupiter Distance (AU): ";
  339.   color 5,0
  340.   print rr;"    "
  341.  
  342.   '..headings
  343.   locate 7,startcol
  344.   color 7,0
  345.   print " Io"
  346.   locate 7,startcol+15
  347.   color 6,0
  348.   print " Europa"
  349.   locate 7,startcol+30
  350.   color 7,0
  351.   print " Ganymede"
  352.   locate 7,startcol+45
  353.   color 6,0
  354.   print " Callisto"
  355.   print
  356.  
  357.   '..satellite's X coordinate
  358.   locate csrlin,10:print "X: ";
  359.   col=startcol
  360.   for moon=1 to 4
  361.     locate csrlin,col
  362.     color moon
  363.     print x(moon);"    ";
  364.     col=col+15 
  365.   next
  366.  
  367.   '..satellite's Y coordinate
  368.   print
  369.   locate csrlin,10:print "Y: ";
  370.   col=startcol
  371.   for moon=1 to 4
  372.     locate csrlin,col
  373.     color moon
  374.     print y(moon);"    ";
  375.     col=col+15  
  376.   next
  377. END SUB
  378.  
  379. '..main
  380. if arg$(1)="?" then 
  381.   print
  382.   print "usage: Jovian date start duration interval"
  383.   print
  384.   print "       where: date is of the form dd mm yyyy"
  385.   print
  386.   print "              start, duration and interval"
  387.   print "              consist of hh mm ss "
  388.   print
  389.   print "              Date and Time are taken to be UT.
  390.   print        
  391.   stop
  392. end if
  393.  
  394. if argcount=12 then  
  395.   dd = val(arg$(1))
  396.   mm = val(arg$(2))
  397.   yy = val(arg$(3))
  398.  
  399.   hrs = val(arg$(4))
  400.   mins = val(arg$(5))
  401.   secs = val(arg$(6))
  402.  
  403.   hr_dur = val(arg$(7))
  404.   min_dur = val(arg$(8))
  405.   sec_dur = val(arg$(9))
  406.  
  407.   hr_int = val(arg$(10))
  408.   min_int = val(arg$(11))
  409.   sec_int = val(arg$(12))
  410. else
  411.   window 1,"Jovian Satellite Simulation",(0,0)-(640,70),0
  412.   font "topaz",8
  413.   print
  414.   print "Start Date (UT):"
  415.   print
  416.   input "Enter day   ",dd
  417.   input "Enter month ",mm
  418.   input "Enter year  ",yy
  419.  
  420.   cls
  421.   print
  422.   print "Start Time (UT):"
  423.   print
  424.   input "Enter hours   ",hrs
  425.   input "Enter minutes ",mins
  426.   input "Enter seconds ",secs
  427.  
  428.   cls
  429.   print
  430.   print "Duration:"
  431.   print
  432.   input "Enter hours   ",hr_dur
  433.   input "Enter minutes ",min_dur
  434.   input "Enter seconds ",sec_dur
  435.  
  436.   cls
  437.   print
  438.   print "Interval:"
  439.   print
  440.   input "Enter hours   ",hr_int
  441.   input "Enter minutes ",min_int
  442.   input "Enter seconds ",sec_int
  443.   window close 1
  444.  
  445.   '..An interval of 6 minutes is the finest resolution 
  446.   '..possible with single-precision floating-point math!
  447.   if hr_int=0 and min_int<6 then min_int=6
  448. end if
  449.  
  450. day_fraction = decimal_hours(hrs,mins,secs) / 24!
  451. JDE = JulianDay(dd,mm,yy) + day_fraction
  452.  
  453. duration = decimal_hours(hr_dur,min_dur,sec_dur) / 24!
  454. end_point = JDE+duration
  455.  
  456. interval = decimal_hours(hr_int,min_int,sec_int) / 24!
  457.  
  458. screen 1,640,200,4,2
  459. window 1,"Jovian Satellite Simulation",(0,0)-(640,200),0,1
  460.  
  461. font "topaz",8
  462.  
  463. palette 0,0,0,0        '..black
  464. palette 1,1,.73,0    '..orange (Io)
  465. palette 2,1,.87,.73    '..tan (Europa, Jupiter)
  466. palette 3,.93,.2,0    '..fire engine red (Ganymede) 
  467. palette 4,.8,.6,.53    '..brown (Callisto) 
  468. palette 5,1,1,1        '..white
  469. palette 6,0,.93,.87    '..aqua
  470. palette 7,.33,.87,0    '..green
  471. palette 8,0,0,1        '..blue
  472. palette 9,1,1,.13    '..yellow
  473.  
  474. MENU 1,0,1,"Project"
  475. MENU 1,1,1,"About..."
  476. MENU 1,2,1,"Quit","Q"
  477.  
  478. ON MENU GOSUB handle_menu
  479. MENU ON
  480.  
  481. '..N,S,E,W markers
  482. color 9,0
  483. locate 15,5
  484. print "E"
  485. locate 15,76
  486. print "W"
  487. locate 12,41
  488. print "N"
  489. locate 18,41
  490. print "S"
  491.  
  492. '..initialise lastX array
  493. for moon=1 to 4
  494.   lastx(moon) = -99
  495. next
  496.  
  497. first=true
  498.  
  499. repeat
  500.   JovianEphemeris(JDE)
  501.   AngleFromInfConj
  502.   DistFromCenter
  503.  
  504.   if not first then 
  505.     RemoveMoons
  506.   else 
  507.     first=false
  508.   end if
  509.  
  510.   for moon=1 to 4
  511.     if not first then lastx(moon) = x(moon)
  512.     calc_x_y(moon)
  513.   next
  514.  
  515.   ShowJovianSpace
  516.  
  517.   DisplayData
  518.  
  519.   JDE = JDE + interval
  520. until JDE > end_point
  521.  
  522. '..await menu selection
  523. while -1:sleep:wend
  524.  
  525. '..trap handlers
  526. handle_menu:
  527.   the_menu = MENU(0)
  528.   the_item = MENU(1)
  529.  
  530.   if the_menu = 1 and the_item = 1 then
  531.     result = MsgBox("Written in ACE by David Benn","OK")
  532.   end if
  533.  
  534.   if the_menu = 1 and the_item = 2 then quit
  535. RETURN
  536.  
  537. quit:
  538.   window close 1
  539.   screen close 1
  540. END
  541.